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The probability density function (PDF) of a global measure in a large class of highly correlated 
systems has been suggested to be of the same functional form. Here, we identify the analytical form 
of the PDF of one such measure, the order parameter in the low temperature phase of the 2D-XY 
model. We demonstrate that this function describes the fluctuations of global quantities in other 
correlated, equilibrium and non-equilibrium systems. These include a coupled rotor model, Ising 
and percolation models, models of forest fires, sand-piles, avalanches and granular media in a self 
organized critical state. We discuss the relationship with both Gaussian and extremal statistics. 
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Self similarity is an important feature of the natural 
world. It arises in strongly correlated many body sys- 
tems when fluctuations over all scales from a microscopic 
length a to a diverging correlation length £ lead to the ap- 
pearence of "anomalous dimension" jl| and fractal prop- 
erties. However, even in an ideal world the divergence of 
of £ must ultimately be cut off by a macroscopic length 
L, allowing the definition of a range of scales between a 
and L, over which the anomalous behaviour can occur. 
Such systems are found, for example, in critical phenom- 
ena, m Self-Organized Criticality |.| or in turbulent 
flow problems. By analogy with fluid mechanics we shall 
call these finite size critical systems "inertial systems" 
and the range of scales between a and L the "inertial 
range" . One of the anomalous statistical properties of in- 
ertial systems is that, whatever their size, they can never 
be divided into mesoscopic regions that are statistically 
independent. As a result they do not satisfy the basic 
criterion of the central limit theorem and one should not 
necessarily expect global, or spatially averaged quanti- 
ties to have Gaussian fluctuations about the mean value. 
In Ref. Q(BHP) it was demonstrated that two of these 
systems, a model of finite size critical behaviour and a 
steady state in a closed turbulent flow experiment, share 
the same non-Gaussian probability distribution function 
(PDF) for fluctuations of global quantities. Consequently 
it was proposed that these two systems - so utterly dis- 
similar in regards to their microscopic details - share the 
same statistics simply because they are critical. If this 
is the case, one should then be able to describe turbu- 
lence as a finite-size critical phenomenon, with an effec- 
tive "universality class". As, however, turbulence and 
the magnetic model arc very unlikely to share the same 
universality class, it was implied that the differences that 
separate critical phenomena into universality classes rep- 
resent at most a minor perturbation on the functional 
form of the PDF. In this paper, to test this proposition, 



we determine the functional form of the BHP fluctuation 
spectrum and show that it indeed applies to a large class 
of inertial systems H) ■ 

The magnetic model studied by BHP, the spin wave 
limit to the two dimensional XY (2D-XY) model, is de- 
fined by the harmonic Hamiltonian 
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where J is the near neighbour exchange constant for an- 
gular variables Oi that occupy a square lattice with peri- 
odic boundary conditions. The magnetization is defined 
as to — l/Nj^i cos(#j — &), where 9 is the instantaneous 
mean orientation. This model is critical at all tempera- 
tures and for an infinite system has algebraic correlations 
on all length scales. In the finite system the lattice con- 
stant a and the system sizes L = ayN define a natural 
inertial range. The model can be diagonalized in Fourier 
space, which makes it very convenient for analytical work. 
The PDF of the magnetization P(m) can be expressed 
as the Fourier transform of a sum over its moments. In 
Ref. 0] it was shown that the moments are given by 
Mn = 9n (32/2) _n / 2 fj™, where a 1 is the variance and the 
gk(k = 2,3,4...) are sums related to the lattice Green 
function in Fourier space G(q): gk — J2q G(q) k /N k . 
The fact that //„ oc /i™ means that a change of N or 
T is equivalent to a linear transformation of the variate 
to; hence, the PDF can be expressed in a universal form. 
As shown in || the moment series can be resummed to 
give the following expression, exact to leading order in 
N: 
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Here y = (m ~ (to) ) / a and (to) is the mean of the dis- 
tribution. Including only 52 in (Q) would give a Gaus- 
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sian PDF with variance a 1 . However the terms for 
k > 2 cannot be neglected and II (y) = crP(y) is a non- 
Gaussian, universal function, independent of both the 
size of the system and the temperature. Without loss 
of generality one can make the quadratic approximation 
m =1 — ~ 0) 2 /2N which allows us to transform 

Eqn. (j^) to a form suitable for numerical integration [fij : 
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>.'7/-\/y - i^TrG/iV — ^Trlog(l - ixG/iV) 



(Here the trace Tr of any function of G is defined as the 
sum for q ^ of the same function of G(q), which can 
be simply proved by making a Taylor expansion of the 
function near G = 0.) In order to make an accurate test 
of this expression we have performed a high resolution 
molecular dynamics simulation of P(m). Fig. 1 compares 
the integrated equation || with data for a system of 1024 
classical rotors integrated over 10 8 molecular dynamics 
time steps in the low temperature phase. The agreement 
is globally excellent, particularly in the wings of the dis- 
tribution and along the exponential tail for fluctuations 
below the mean. 

The asymptotic values of H(y) are related to the saddle 
points of the integrand in (||). We find 



U{y) oc 1 2/1 exp(--by); for y < 1 



Tl(y) cx exp(-|e b( ^ s) ); for y > 1 
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where b = ^^fg%f2 1.105 and s = 0.745. These forms 
give the correct asymptotic gradients of the molecular 
dynamics data on logarithmic and double logarithmic 
scales. The asymptotic forms are an accurate approxima- 
tion to Eqn. (Q) for large |y|; however deviations from the 
asymptotes are important over most of the physical range 
of y, which is typically limited to log|y| ~ O(l). Equa- 
tions (|]) and (||) serve as a guide to finding a good ap- 
proximation to the functional form of Tl(y) in this range. 
To do this, we observe that the factor of y in (^) can be 
regarded as constant in this regime, which along with (^) 
and (@) immediately suggests the form 



H(y) = K 
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This function must obey the three conditions of unit area, 
zero mean and unit variance, which fixes 6, s and K 
to values slightly different to those found analytically: 
b = 0.938, s = 0.374, K = 2.14. An alternative approach 
is to choose the parameters in the generalized function 



jy e a(b(y-s)-e b< - v s) ) guc j 1 fa&t the first four Fourier coef- 
ficients match Eqn.(||). In this case we find a = 1.58, 
K = 2.16, b = 0.934, s — 0.373, in satisfying agreement 
with the previous estimates. The ratios of the higher or- 
der Fourier coefficients differ from unity only very slowly, 
showing that Eqn. (||) is an accurate approximation to 
n(y). This is directly confirmed by plotting Eqn. (^) ver- 
sus the molecular dynamics and exact results in Fig. 1, 
where the fit is seen to be of extremely high precision. 
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FIG. 1. The PDF II(y) as found from a numerical Fourier 
transform of Eqn. (^|) (long dash), from Eqn. ([]) (solid) and 
by Molecular Dynamics simulation of a system of N = 1024 
classical rotors (dotted). 

We now test the idea that the BHP fluctuation spec- 
trum of the form of Eqn. ^ is exhibited by many types of 
inertial system. Fig. 2 shows the numerically simulated 
PDF of global quantities in several equilibrium and non- 
equilibrium models. The equilibrium models include the 
2D Ising model at a temperature T*(N) just below the 
critical temperature and a 2D site percolation model on 
a square lattice for a site occupation probability P*(N) 
just above the percolation transition. The numerical re- 
sults refer to the fluctuations of the absolute value of 
the magnetization and the fluctuations in the size of the 
spanning cluster respectively. The non-equilibrium mod- 
els are of the type that when driven slowly enter a scale 
free or critical steady state defined as Self-Organized 
Criticality (SOC) ^|J3|. Here the global quantity is es- 
sentially a dissipation rate that fluctuates about a well 
defined mean value in the steady state. Details of the 
individual SOC models are as follows. (1) The auto ig- 
niting forest fire model ^ consists of "trees" planted at 
random on the vacant sites of a square lattice with prob- 
ability p. In each time step the age Ti of a tree on site 
i is incremented by one unit. When Ti = T max the tree 
ignites and Ti is reset to zero. Trees can also catch fire 
by being nearest neighbour to a site on fire. The energy, 
or wood stored in a tree, is proportional to T, and the 
figure shows the PDF of the total energy dissipated in 
fires at each time step. (2) In the Bak-Tang-Wiesenfeld 
(BTW) sandpile model M a dynamical variable Ei is de- 
fined on lattice site i. The model is driven by adding 
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units of the Enfield to randomly selected sites. When 
Ei > E max the site variable is decreased by E max and 
the ^-variable of the z neighbour sites are increased by 
E m ax/ z. One or more of the neighbour sites may then 
aquirc an E- value larger than E max and an avalanche is 
induced. The PDF shown refers to the fluctuations in the 
instantaneous number of relaxing sites. (3) In the Snep- 
pen depinning model JTo| an interface moves through a 
static random field of pinning forces. The site along the 
interface that experiences the smallest pinning force is 
moved one unit ahead. If the local slope Sj exceeds 1, 
then the neighbouring sites are moved one unit ahead 
until all Si < 1. We call such a sequence of updates 
a micro-avalanche and calculate the PDF of the sum of 
areas covered by the progressing interface during an in- 
tegral time scale T dependent on system size. (4) The 
model for granular media is a 'Tetris-like' 2D lattice gas 
ensemble of anisotropic particles settling under gravity 
in a finite box Hj. Due to the geometrical frustration, 
the total mass varies from one realization of the filling 
process to another. The PDF for fluctuations in bulk 
density of the particles is shown. 
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FIG. 2. Fluctuation spectra in equilibrium systems (lower 
curves): Ising (o), percolation (A). The central curve (V) 
correspond to the fluctuation spectrum of the correlated ex- 
tremal process - see text. The upper curves are the PDFs for 
the auto igniting forest fire model (x), the Sneppen depin- 
ning model (+), the granular media model (□) and the BTW 
sandpile model (*). The dotted line corresponds to Eqn. (|^). 
For clarity, the sets of curves are shifted downwards by 1.5 in 
log units. 

Refering to Fig. 2, the data sets for all models fall 
close to the BHP form, Eqn. (||). In the equilibrium 
models (lower curves) the self similarity is expected at 
the system-size dependent critical temperature T* (N) or 
percolation probability P*(N) only. The PDF for the 2D 
Ising model, for example, is temperature dependent, but 
makes a close approach to the BHP form around T* (N) . 
We believe that the remaining deviations for fluctuations 
above the mean are due to the limited inertial range for 
the system sizes studied. For the non-equilibrium sys- 
tems (upper curves) some of the data sets also show some 



deviation. This may simply be due to poor statistics, as 
the deviations on the two sides of the mean are related 
by the constraints of normalisation. In this respect, we 
note that extremely good statistics were required to get 
a really satisfactory fit to Eqn. (^|) for the 2D-XY model, 
while for more limited data sets systematic deviations 
occurred. 

Nevertheless it is clear that, to leading order, Eqn. (|^) 
correctly gives the behaviour of the global fluctuation 
spectrum in all these systems, independently of the de- 
tails of the each example. We propose that this is a con- 
sequence of the systems sharing the properties of finite 
size, strong correlations and self similarity. 

To clarify this proposition, we return to our calcula- 
tion on the 2D-XY model. One can see explicitly that 
the BHP spectrum occurs through appearance of anoma- 
lous dimension and contributions on all length scales of 
the intertial range. The magnetization can be written, 
within the quadratic approximation, m — 1 — E q rn q , 
where the m q s are the amplitudes from the individual 
spin wave modes. These are statistically independent 
positive variates with PDF 



P(m q ) = ^^L_m- 1 / 2 exp(-(3Nq 2 m q ), (7) 

whose mean and standard deviation therefore scale with 
q~ 2 . The "softest" modes have wave vector q = 2n/L 
and hence, by themselves, make contributions of 0(1) 
to m, while the modes on the zone boundary with 
q = it /a have only microscopic amplitude. The mo- 
ments of P(m) are determined by the mean magneti- 
zation. This is proportional to the integral over all con- 
tributions: ~ j^J L q~ 2 n(q)dq where n(q) ~ q^ 1 is the 
density of states. In one dimension the integral depends 
only on the lower limit 2n/L and only the soft modes 
count, while in three dimensions only the upper limit 
it /a is important and the multitude of modes near the 
zone boundary dominate the sum. In two-dimensions 
however, both limits of the integral are required and a 
detailed calculation gives (m) = 1 — 77/2 log (CL/a), with 
C = 1.87 @ and critical exponent 77 = T/2ttJ. The rel- 
evance of fluctuations over all length scales of the zone 
therefore leads to the "anomalous" term log {L/a) and it 
ensures that the system cannot be cut into statistically 
independent parts. 

The spinwave approximation to the XY-model is ex- 
actly equivalent to the Edwards- Wilkinson (EW) model 
of a growing interface in steady state |p^| , with the square 
of the interface width w equal to the sum over the ampli- 
tudes m q : w = X)g m 9- The fluctuations in the width 
of the interface has been studied by Foltin et al. |l3| for 
the ID case and by Racz and Plischke (lj] for the 2D in- 
terface. The BHP spectrum is found for the critical two 
dimensional case only and our calculation can be consid- 
ered the completion of the study in Ref. . 
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The functional form of Eq n. (p| ) suggests a relationship 
to Gumbel's first asymptote |15| for extreme value statis- 
tics, which have recently been discussed in relation to 
turbulence in one dimension jn| . The form but with 
a taking integer values, where a = 1,2,3... would corre- 
spond to the PDF for the first, second, third.... largest of 
the N random numbers. However, the exponent a = n/2 
suggests, as we have argued, that the fluctuations in 
m are not dominated by single independent variables. 
Rather, the analytic derivation of Eqn. (^) shows that if 
extreme value statistics are involved they must be related 
to the statistics of some emergent coherent collective ex- 
citation of the system. This is born out in the simulations 
of the Ising model and of all the SOC models studied. In 
the Ising model, it is found that both the full magnetiza- 
tion and the contribution to the magnetisation from the 
largest connected cluster of parallel spins give the same 
PDF, within numerical error. For the Sneppen model the 
PDF of the sum over avalanches and that of the largest 
avalanche during the integral time T, are both found to 
be of the BHP form, even though these quantities are 
not related by a simple scale. If the avalanches appear- 
ing during time T were uncorrelated one would expect 
that the PDF for the largest avalanche would be Gum- 
bel's asymptote with a = 1. The modification towards 
our form indicates therefore that there are correlations 
between events during the period T. To test this idea, 
we have studied the PDF of the extreme values taken 
from sets of linearly correlated variables. The process 
consists of generating a vector \ — (Xii ■■■Xn) where \i 
are all independent and exponentially distributed. The 
maximum signal is obtained as ^ m ax — rnax{£i, £/v} 
where the vector £ = and M is an TV x N matrix with 
random but fixed elements. The resulting spectrum, the 
central data set in Fig. 2, is found to be very close to the 
BHP spectrum. 

In conclusion, our results infer that the non-Gaussian 
PDF of a global quantity in a critical system is a con- 
sequence of finite size, strong correlations and self simi- 
larity, and is independent of universality class to leading 
order. Clearly many more studies of this point are re- 
quired. Non-linearity does not appear to be an essential 
feature, over and above the necessity, in a closed sys- 
tem, to couple the elementary degrees of freedom. In- 
deed, if non-linearity were essential, it would seem im- 
possible that the linear spin-wave theory could capture 
the fluctuations in the turbulence experiment Q]. Racz 
and Plischke Jl4|] have studied a series of linear and non- 
linear models for growing interfaces. All show anisotropic 
PDFs for the interface width, with long tails in qualita- 
tive agreement with our data. It would very interesting 
to examine these models in detail to see how, in a con- 
trolled environment the strength of the non-linearity af- 
fects the form of the fluctuations. Finally, it seems that a 
relationship exists between the BHP curve and extremal 



statistics. Although we have shown that the BHP be- 
havior is not simply due to extreme values of the sta- 
tistically independent degrees of freedom of the 2D-XY 
model, extreme values do appear to dominate the real 
space coherent structures that are excited in the critical 
(Ising model) or self-organized critical (Sneppen model) 
state. 

Our findings thus establish a completely new and gen- 
eral consequence of self similarity and they open the door 
to numerous studies that could lead to a unified global 
description of aspects of equilibrium and non-equilibrium 
behaviour. 
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